{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m    Updating\u001b[22m\u001b[39m registry at `C:\\Users\\danil\\.julia\\registries\\General`\n",
      "\u001b[32m\u001b[1m    Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Manifest.toml`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Manifest.toml`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Manifest.toml`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Manifest.toml`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `C:\\Users\\danil\\.julia\\environments\\v1.8\\Manifest.toml`\n"
     ]
    }
   ],
   "source": [
    "import Pkg\n",
    "Pkg.add(\"Plots\")\n",
    "Pkg.add(\"Distributions\")\n",
    "Pkg.add(\"Statistics\")\n",
    "Pkg.add(\"StatsBase\")\n",
    "Pkg.add(\"Random\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "using Distributions\n",
    "using Statistics\n",
    "using StatsBase\n",
    "using Random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TaskLocalRNG()"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setting the random seed\n",
    "Random.seed!(19013)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100×1 Matrix{Float64}:\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " ⋮\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simulations of one-period expected posterior as a function of q\n",
    "\n",
    "gamma = 0.8\n",
    "gammahat = 0.5\n",
    "\n",
    "p1 = 0.05\n",
    "p2 = 0.1\n",
    "p3 = 0.2\n",
    "p4 = 0.4\n",
    "\n",
    "n = 1\n",
    "dL = 10\n",
    "dR = 10\n",
    "\n",
    "qvec = range(0.505,length=100,stop=0.995)\n",
    "\n",
    "muR1 = zeros(100,1)\n",
    "muR2 = zeros(100,1)\n",
    "muR3 = zeros(100,1)\n",
    "muR4 = zeros(100,1)\n",
    "\n",
    "muL1 = zeros(100,1)\n",
    "muL2 = zeros(100,1)\n",
    "muL3 = zeros(100,1)\n",
    "muL4 = zeros(100,1)\n",
    "\n",
    "mu1 = zeros(100,1)\n",
    "mu2 = zeros(100,1)\n",
    "mu3 = zeros(100,1)\n",
    "mu4 = zeros(100,1)\n",
    "\n",
    "pratio1 = zeros(100,1)\n",
    "pratio2 = zeros(100,1)\n",
    "pratio3 = zeros(100,1)\n",
    "pratio4 = zeros(100,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\shortrun_posterior_n1a10b10.png\""
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for k in 1:1:100\n",
    "    for i in 0:1:n\n",
    "        for j in 0:1:(dL+dR)\n",
    "            muR1[k] = muR1[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p1/(p1+(1-p1)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR2[k] = muR2[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p2/(p2+(1-p2)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR3[k] = muR3[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p3/(p3+(1-p3)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR4[k] = muR4[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p4/(p4+(1-p4)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    for i in 0:1:n\n",
    "        for j in 0:1:(dL+dR)\n",
    "            muL1[k] = muL1[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p1/(p1+(1-p1)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL2[k] = muL2[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p2/(p2+(1-p2)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL3[k] = muL3[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p3/(p3+(1-p3)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL4[k] = muL4[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p4/(p4+(1-p4)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    mu1[k] = p1*muL1[k]+(1-p1)*muR1[k]\n",
    "    mu2[k] = p2*muL2[k]+(1-p2)*muR2[k]\n",
    "    mu3[k] = p3*muL3[k]+(1-p3)*muR3[k]\n",
    "    mu4[k] = p4*muL4[k]+(1-p4)*muR4[k]\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    pratio1[k] = mu1[k]/p1\n",
    "    pratio2[k] = mu2[k]/p2\n",
    "    pratio3[k] = mu3[k]/p3\n",
    "    pratio4[k] = mu4[k]/p4\n",
    "end\n",
    "\n",
    "onelevel = ones(100,1)\n",
    "\n",
    "plot(qvec,[pratio1 pratio2 pratio3 pratio4 onelevel], label=[\"π=0.05\" \"π=0.1\" \"π=0.2\" \"π=0.4\" \"\"], xlabel = \"q\", ylabel = \"E[μ(s₁)] / π\", legend=:bottomright, legendfontsize=12, labelfontsize=14, tickfontsize=12, xlimits=[0.5,1], ylimits=[0.88,1.06], linecolor=[:blue :orange :green :red :black], linestyle=[:solid :dash :dot :dashdot :dash], size=(800,600), dpi=600)\n",
    "savefig(\"shortrun_posterior_n1a10b10\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100×1 Matrix{Float64}:\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " ⋮\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma = 0.8\n",
    "gammahat = 0.5\n",
    "\n",
    "p1 = 0.05\n",
    "p2 = 0.1\n",
    "p3 = 0.2\n",
    "p4 = 0.4\n",
    "\n",
    "n = 1\n",
    "dL = 11\n",
    "dR = 10\n",
    "\n",
    "qvec = range(0.505,length=100,stop=0.995)\n",
    "\n",
    "muR1 = zeros(100,1)\n",
    "muR2 = zeros(100,1)\n",
    "muR3 = zeros(100,1)\n",
    "muR4 = zeros(100,1)\n",
    "\n",
    "muL1 = zeros(100,1)\n",
    "muL2 = zeros(100,1)\n",
    "muL3 = zeros(100,1)\n",
    "muL4 = zeros(100,1)\n",
    "\n",
    "mu1 = zeros(100,1)\n",
    "mu2 = zeros(100,1)\n",
    "mu3 = zeros(100,1)\n",
    "mu4 = zeros(100,1)\n",
    "\n",
    "pratio1 = zeros(100,1)\n",
    "pratio2 = zeros(100,1)\n",
    "pratio3 = zeros(100,1)\n",
    "pratio4 = zeros(100,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\shortrun_posterior_n1a11b10.png\""
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for k in 1:1:100\n",
    "    for i in 0:1:n\n",
    "        for j in 0:1:(dL+dR)\n",
    "            muR1[k] = muR1[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p1/(p1+(1-p1)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR2[k] = muR2[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p2/(p2+(1-p2)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR3[k] = muR3[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p3/(p3+(1-p3)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR4[k] = muR4[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p4/(p4+(1-p4)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    for i in 0:1:n\n",
    "        for j in 0:1:(dL+dR)\n",
    "            muL1[k] = muL1[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p1/(p1+(1-p1)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL2[k] = muL2[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p2/(p2+(1-p2)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL3[k] = muL3[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p3/(p3+(1-p3)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL4[k] = muL4[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p4/(p4+(1-p4)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    mu1[k] = p1*muL1[k]+(1-p1)*muR1[k]\n",
    "    mu2[k] = p2*muL2[k]+(1-p2)*muR2[k]\n",
    "    mu3[k] = p3*muL3[k]+(1-p3)*muR3[k]\n",
    "    mu4[k] = p4*muL4[k]+(1-p4)*muR4[k]\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    pratio1[k] = mu1[k]/p1\n",
    "    pratio2[k] = mu2[k]/p2\n",
    "    pratio3[k] = mu3[k]/p3\n",
    "    pratio4[k] = mu4[k]/p4\n",
    "end\n",
    "\n",
    "onelevel = ones(100,1)\n",
    "\n",
    "plot(qvec,[pratio1 pratio2 pratio3 pratio4 onelevel], label=[\"π=0.05\" \"π=0.1\" \"π=0.2\" \"π=0.4\" \"\"], xlabel = \"q\", ylabel = \"E[μ(s₁)] / π\", legend=:bottomright, legendfontsize=12, labelfontsize=14, tickfontsize=12, xlimits=[0.5,1], ylimits=[0.88,1.06], linecolor=[:blue :orange :green :red :black], linestyle=[:solid :dash :dot :dashdot :dash], size=(800,600), dpi=600)\n",
    "savefig(\"shortrun_posterior_n1a11b10\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100×1 Matrix{Float64}:\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " ⋮\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma = 0.8\n",
    "gammahat = 0.5\n",
    "\n",
    "p1 = 0.05\n",
    "p2 = 0.1\n",
    "p3 = 0.2\n",
    "p4 = 0.4\n",
    "\n",
    "n = 1\n",
    "dL = 12\n",
    "dR = 10\n",
    "\n",
    "qvec = range(0.505,length=100,stop=0.995)\n",
    "\n",
    "muR1 = zeros(100,1)\n",
    "muR2 = zeros(100,1)\n",
    "muR3 = zeros(100,1)\n",
    "muR4 = zeros(100,1)\n",
    "\n",
    "muL1 = zeros(100,1)\n",
    "muL2 = zeros(100,1)\n",
    "muL3 = zeros(100,1)\n",
    "muL4 = zeros(100,1)\n",
    "\n",
    "mu1 = zeros(100,1)\n",
    "mu2 = zeros(100,1)\n",
    "mu3 = zeros(100,1)\n",
    "mu4 = zeros(100,1)\n",
    "\n",
    "pratio1 = zeros(100,1)\n",
    "pratio2 = zeros(100,1)\n",
    "pratio3 = zeros(100,1)\n",
    "pratio4 = zeros(100,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\shortrun_posterior_n1a12b10.png\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for k in 1:1:100\n",
    "    for i in 0:1:n\n",
    "        for j in 0:1:(dL+dR)\n",
    "            muR1[k] = muR1[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p1/(p1+(1-p1)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR2[k] = muR2[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p2/(p2+(1-p2)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR3[k] = muR3[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p3/(p3+(1-p3)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "            muR4[k] = muR4[k] + binomial(n,i)*(1-qvec[k])^i*(qvec[k])^(n-i)*binomial(dL+dR,j)*(1-qvec[k])^j*(qvec[k])^(dL+dR-j)*p4/(p4+(1-p4)*((1-qvec[k])/qvec[k])^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*qvec[k]+(1-gamma)))^(j-dL))\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    for i in 0:1:n\n",
    "        for j in 0:1:(dL+dR)\n",
    "            muL1[k] = muL1[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p1/(p1+(1-p1)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL2[k] = muL2[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p2/(p2+(1-p2)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL3[k] = muL3[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p3/(p3+(1-p3)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "            muL4[k] = muL4[k] + binomial(n,i)*(qvec[k])^i*(1-qvec[k])^(n-i)*binomial(dL+dR,j)*(qvec[k])^j*(1-qvec[k])^(dL+dR-j)*p4/(p4+(1-p4)*((1-qvec[k])/(qvec[k]))^(2*i+j-dR-n)*((gamma*(1-qvec[k])+(1-gamma))/(gamma*(qvec[k])+(1-gamma)))^(j-dL))\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    mu1[k] = p1*muL1[k]+(1-p1)*muR1[k]\n",
    "    mu2[k] = p2*muL2[k]+(1-p2)*muR2[k]\n",
    "    mu3[k] = p3*muL3[k]+(1-p3)*muR3[k]\n",
    "    mu4[k] = p4*muL4[k]+(1-p4)*muR4[k]\n",
    "end\n",
    "\n",
    "for k in 1:1:100\n",
    "    pratio1[k] = mu1[k]/p1\n",
    "    pratio2[k] = mu2[k]/p2\n",
    "    pratio3[k] = mu3[k]/p3\n",
    "    pratio4[k] = mu4[k]/p4\n",
    "end\n",
    "\n",
    "onelevel = ones(100,1)\n",
    "\n",
    "plot(qvec,[pratio1 pratio2 pratio3 pratio4 onelevel], label=[\"π=0.05\" \"π=0.1\" \"π=0.2\" \"π=0.4\" \"\"], xlabel = \"q\", ylabel = \"E[μ(s₁)] / π\", legend=:bottomright, legendfontsize=12, labelfontsize=14, tickfontsize=12, xlimits=[0.5,1], ylimits=[0.88,1.06], linecolor=[:blue :orange :green :red :black], linestyle=[:solid :dash :dot :dashdot :dash], size=(800,600), dpi=600)\n",
    "savefig(\"shortrun_posterior_n1a12b10\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simulations of beliefs' evolution over many periods\n",
    "\n",
    "# Parameters and variables\n",
    "\n",
    "gamma = 0.95\n",
    "gammahat = 0.9*gamma\n",
    "p = 0.5\n",
    "\n",
    "# Friend networks\n",
    "n = 5\n",
    "dMaj = 6\n",
    "dMin = 4\n",
    "\n",
    "T = 300\n",
    "\n",
    "sMaj1 = 0\n",
    "sMin1 = 0\n",
    "sn1 = 0\n",
    "sAmaj1 = 0\n",
    "sAmin1 = 0\n",
    "sAn1 = 0\n",
    "\n",
    "sMaj2 = 0\n",
    "sMin2 = 0\n",
    "sn2 = 0\n",
    "sAmaj2 = 0\n",
    "sAmin2 = 0\n",
    "sAn2 = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.505984984984985"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculation of q_LR (variable qbar in this code)\n",
    "\n",
    "qvec = range(0.501,length=1000,stop=0.999)\n",
    "z = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    z[k] = log((gammahat*(1-qvec[k])+(1-gammahat))/(gammahat*qvec[k]+(1-gammahat)))/log((1-qvec[k])/qvec[k])\n",
    "end\n",
    "qb = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    qb[k] = abs(1/2+((gamma-(2-gamma)*z[k])*(dMaj-dMin))/(2*gamma*(2*(1+n)+(dMaj+dMin)*(1+z[k])))-qvec[k])\n",
    "end\n",
    "diff = 0\n",
    "ind = 0\n",
    "diff, ind = findmin(qb)\n",
    "qbar = qvec[ind]\n",
    "\n",
    "qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5017954954954955"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the low-q case\n",
    "\n",
    "q = 0.7*0.5+(1-0.7)*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beliefs with misperception\n",
    "\n",
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beliefs without misperception\n",
    "\n",
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\5n6a4b_lowq_nomisp.png\""
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: small friend network (15), small imbalance (6:4), low q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-6a-4b, misperception, low q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n6a4b_lowq_misp\")\n",
    "\n",
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-6a-4b, no misperception, low q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n6a4b_lowq_nomisp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5041894894894895"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the high-q case\n",
    "\n",
    "q = 0.3*0.5+(1-0.3)*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beliefs with misperception\n",
    "\n",
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beliefs without misperception\n",
    "\n",
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\5n6a4b_highq_nomisp.png\""
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: small friend network (15), small imbalance (6:4), high q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-6a-4b, misperception, high q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n6a4b_highq_misp\")\n",
    "\n",
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-6a-4b, no misperception, high q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n6a4b_highq_nomisp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Parameters and variables\n",
    "\n",
    "gamma = 0.95\n",
    "gammahat = 0.9*gamma\n",
    "p = 0.5\n",
    "\n",
    "# Friend networks\n",
    "\n",
    "n = 5\n",
    "dMaj = 9\n",
    "dMin = 1\n",
    "\n",
    "T = 300\n",
    "\n",
    "sMaj1 = 0\n",
    "sMin1 = 0\n",
    "sn1 = 0\n",
    "sAmaj1 = 0\n",
    "sAmin1 = 0\n",
    "sAn1 = 0\n",
    "\n",
    "sMaj2 = 0\n",
    "sMin2 = 0\n",
    "sn2 = 0\n",
    "sAmaj2 = 0\n",
    "sAmin2 = 0\n",
    "sAn2 = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5239309309309309"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculation of q_LR (variable qbar in this code)\n",
    "\n",
    "qvec = range(0.501,length=1000,stop=0.999)\n",
    "z = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    z[k] = log((gammahat*(1-qvec[k])+(1-gammahat))/(gammahat*qvec[k]+(1-gammahat)))/log((1-qvec[k])/qvec[k])\n",
    "end\n",
    "qb = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    qb[k] = abs(1/2+((gamma-(2-gamma)*z[k])*(dMaj-dMin))/(2*gamma*(2*(1+n)+(dMaj+dMin)*(1+z[k])))-qvec[k])\n",
    "end\n",
    "diff = 0\n",
    "ind = 0\n",
    "diff, ind = findmin(qb)\n",
    "qbar = qvec[ind]\n",
    "\n",
    "qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5071792792792793"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the low-q case\n",
    "\n",
    "q = 0.7*0.5 + (1-0.7)*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\5n9a1b_lowq_nomisp.png\""
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: small echo chamber (15), large imbalance (9:1), low q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-9a-1b, misperception, low q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n9a1b_lowq_misp\")\n",
    "\n",
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-9a-1b, no misperception, low q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n9a1b_lowq_nomisp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5167516516516516"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the high-q case\n",
    "\n",
    "q = 0.3*0.5+(1-0.3)*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\5n9a1b_highq_nomisp.png\""
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: small friend network (15), large imbalance (9:1), high q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-9a-1b, misperception, high q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n9a1b_highq_misp\")\n",
    "\n",
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 5n-9a-1b, no misperception, high q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"5n9a1b_highq_nomisp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Parameters and variables\n",
    "\n",
    "gamma = 0.95\n",
    "gammahat = 0.9*gamma\n",
    "p = 0.5\n",
    "\n",
    "# Friend networks\n",
    "\n",
    "n = 17\n",
    "dMaj = 24\n",
    "dMin = 16\n",
    "\n",
    "T = 300\n",
    "\n",
    "sMaj1 = 0\n",
    "sMin1 = 0\n",
    "sn1 = 0\n",
    "sAmaj1 = 0\n",
    "sAmin1 = 0\n",
    "sAn1 = 0\n",
    "\n",
    "sMaj2 = 0\n",
    "sMin2 = 0\n",
    "sn2 = 0\n",
    "sAmaj2 = 0\n",
    "sAmin2 = 0\n",
    "sAn2 = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5064834834834835"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculation of q_LR (variable qbar in this code)\n",
    "\n",
    "qvec = range(0.501,length=1000,stop=0.999)\n",
    "z = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    z[k] = log((gammahat*(1-qvec[k])+(1-gammahat))/(gammahat*qvec[k]+(1-gammahat)))/log((1-qvec[k])/qvec[k])\n",
    "end\n",
    "qb = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    qb[k] = abs(1/2+((gamma-(2-gamma)*z[k])*(dMaj-dMin))/(2*gamma*(2*(1+n)+(dMaj+dMin)*(1+z[k])))-qvec[k])\n",
    "end\n",
    "diff = 0\n",
    "ind = 0\n",
    "diff, ind = findmin(qb)\n",
    "qbar = qvec[ind]\n",
    "\n",
    "qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.501945045045045"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the low-q case\n",
    "\n",
    "q = 0.7*0.5+0.3*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\17n24a16b_lowq_misp_legend.png\""
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: large friend network (75), small imbalance (6:4), low q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Mean of Alice's beliefs\" \"10th/90th quantile of Alice's beliefs\" \"\" \"Mean of Bob's beliefs\" \"10th/90th quantile of Bob's beliefs\" \"\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:dash :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), dpi=600, legend=:bottomright, ylabel = \"Belief in state A\", xlabel = \"Period, t\", legendfontsize=12, labelfontsize=14, tickfontsize=12)\n",
    "# title = \"Chamber 17n-24a-16b, misperception, low q\"\n",
    "savefig(\"17n24a16b_lowq_misp_legend\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\17n24a16b_lowq_nomisp_legend.png\""
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Mean of Alice*'s beliefs\" \"10th/90th quantile of Alice*'s beliefs\" \"\" \"Mean of Bob*'s beliefs\" \"10th/90th quantile of Bob*'s beliefs\" \"\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:dash :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), dpi=600, legend=:bottomright, ylabel = \"Belief in state A\", xlabel = \"Period, t\", legendfontsize=12, labelfontsize=14, tickfontsize=12)\n",
    "# title = \"Chamber 17n-24a-16b, no misperception, low q\", \n",
    "savefig(\"17n24a16b_lowq_nomisp_legend\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5045384384384384"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the high-q case\n",
    "\n",
    "q = 0.3*0.5+0.7*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\17n24a16b_highq_misp_legend.png\""
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: large friend network (75), small imbalance (6:4), high q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Mean of Alice's beliefs\" \"10th/90th quantile of Alice's beliefs\" \"\" \"Mean of Bob's beliefs\" \"10th/90th quantile of Bob's beliefs\" \"\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:dash :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), dpi=600, legend=:bottomright, ylabel = \"Belief in state A\", xlabel = \"Period, t\", legendfontsize=12, labelfontsize=14, tickfontsize=12)\n",
    "# title = \"Chamber 17n-24a-16b, misperception, high q\", \n",
    "savefig(\"17n24a16b_highq_misp_legend\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\17n24a16b_highq_nomisp_legend.png\""
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Mean of Alice*'s beliefs\" \"10th/90th quantile of Alice*'s beliefs\" \"\" \"Mean of Bob*'s beliefs\" \"10th/90th quantile of Bob*'s beliefs\" \"\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:dash :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), dpi=600, legend=:bottomright, ylabel = \"Belief in state A\", xlabel = \"Period, t\", legendfontsize=12, labelfontsize=14, tickfontsize=12)\n",
    "# title = \"Chamber 17n-24a-16b, no misperception, low q\", \n",
    "savefig(\"17n24a16b_highq_nomisp_legend\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Parameters and variables\n",
    "\n",
    "gamma = 0.95\n",
    "gammahat = 0.9*gamma\n",
    "p = 0.5\n",
    "\n",
    "# Friend networks\n",
    "n = 17\n",
    "dMaj = 36\n",
    "dMin = 4\n",
    "\n",
    "T = 300\n",
    "\n",
    "mu1 = zeros(T+1,1)\n",
    "mu1[1] = p\n",
    "mu2 = zeros(T+1,1)\n",
    "mu2[1] = p\n",
    "\n",
    "sMaj1 = 0\n",
    "sMin1 = 0\n",
    "sn1 = 0\n",
    "sAmaj1 = 0\n",
    "sAmin1 = 0\n",
    "sAn1 = 0\n",
    "\n",
    "sMaj2 = 0\n",
    "sMin2 = 0\n",
    "sn2 = 0\n",
    "sAmaj2 = 0\n",
    "sAmin2 = 0\n",
    "sAn2 = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5264234234234234"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# qbar\n",
    "qvec = range(0.501,length=1000,stop=0.999)\n",
    "z = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    z[k] = log((gammahat*(1-qvec[k])+(1-gammahat))/(gammahat*qvec[k]+(1-gammahat)))/log((1-qvec[k])/qvec[k])\n",
    "end\n",
    "qb = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    qb[k] = abs(1/2+((gamma-(2-gamma)*z[k])*(dMaj-dMin))/(2*gamma*(2*(1+n)+(dMaj+dMin)*(1+z[k])))-qvec[k])\n",
    "end\n",
    "diff = 0\n",
    "ind = 0\n",
    "diff, ind = findmin(qb)\n",
    "qbar = qvec[ind]\n",
    "\n",
    "qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.507927027027027"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the low-q case\n",
    "\n",
    "q = 0.7*0.5+0.3*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\17n36a4b_lowq_nomisp.png\""
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: large friend network (75), large imbalance (9:1), low q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 17n-36a-4b, misperception, low q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"17n36a4b_lowq_misp\")\n",
    "\n",
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 17n-36a-4b, no misperception, low q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"17n36a4b_lowq_nomisp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5184963963963963"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determination of q for the high-q case\n",
    "\n",
    "q = 0.3*0.5+0.7*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1vec_nm = zeros(T+1,10000)\n",
    "mu1vec_nm[1,:] = p*ones(1,10000)\n",
    "mu2vec_nm = zeros(T+1,10000)\n",
    "mu2vec_nm[1,:] = p*ones(1,10000)\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:T\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_nm[t+1,i] = mu1vec_nm[t,i]/(mu1vec_nm[t,i]+(1-mu1vec_nm[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_nm[t+1,i] = mu2vec_nm[t,i]/(mu2vec_nm[t,i]+(1-mu2vec_nm[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gamma*(1-q)+(1-gamma))/(gamma*q+(1-gamma)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\17n36a4b_highq_nomisp.png\""
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Misperception vs no misperception\n",
    "# Parameters: large friend network(75), large imbalance (9:1), high q\n",
    "\n",
    "average_mu1_m = zeros(T+1,1)\n",
    "a10_mu1_m = zeros(T+1,1)\n",
    "a90_mu1_m = zeros(T+1,1)\n",
    "\n",
    "average_mu2_m = zeros(T+1,1)\n",
    "a10_mu2_m = zeros(T+1,1)\n",
    "a90_mu2_m = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_m[t] = mean(mu1vec_m[t,:])\n",
    "    a10_mu1_m[t] = quantile!(mu1vec_m[t,:],0.1)\n",
    "    a90_mu1_m[t] = quantile!(mu1vec_m[t,:],0.9)\n",
    "    \n",
    "    average_mu2_m[t] = mean(mu2vec_m[t,:])\n",
    "    a10_mu2_m[t] = quantile!(mu2vec_m[t,:],0.1)\n",
    "    a90_mu2_m[t] = quantile!(mu2vec_m[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_m a10_mu1_m a90_mu1_m average_mu2_m a10_mu2_m a90_mu2_m]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 17n-36a-4b, misperception, high q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"17n36a4b_highq_misp\")\n",
    "\n",
    "average_mu1_nm = zeros(T+1,1)\n",
    "a10_mu1_nm = zeros(T+1,1)\n",
    "a90_mu1_nm = zeros(T+1,1)\n",
    "\n",
    "average_mu2_nm = zeros(T+1,1)\n",
    "a10_mu2_nm = zeros(T+1,1)\n",
    "a90_mu2_nm = zeros(T+1,1)\n",
    "\n",
    "for t in 1:1:(T+1)\n",
    "    average_mu1_nm[t] = mean(mu1vec_nm[t,:])\n",
    "    a10_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.1)\n",
    "    a90_mu1_nm[t] = quantile!(mu1vec_nm[t,:],0.9)\n",
    "    \n",
    "    average_mu2_nm[t] = mean(mu2vec_nm[t,:])\n",
    "    a10_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.1)\n",
    "    a90_mu2_nm[t] = quantile!(mu2vec_nm[t,:],0.9)\n",
    "end\n",
    "\n",
    "plot_data = [average_mu1_nm a10_mu1_nm a90_mu1_nm average_mu2_nm a10_mu2_nm a90_mu2_nm]\n",
    "plot_labels = [\"Average μᴬ\" \"10% μᴬ\" \"90% μᴬ\" \"Average μᴮ\" \"10% μᴮ\" \"90% μᴮ\"]\n",
    "plot_colors = [:blue :blue :blue :red :red :red]\n",
    "plot_styles = [:solid :dot :dot :solid :dot :dot]\n",
    "plot(plot_data, title = \"Chamber 17n-36a-4b, no misperception, high q\", label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), legend=:outertopright, ylabel = \"Belief in state A\", xlabel = \"Period, t\")\n",
    "savefig(\"17n36a4b_highq_nomisp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5041894894894895"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simulation of climate change opinions (Gallup) with fixed friend networks\n",
    "\n",
    "T = 50\n",
    "n = 5\n",
    "dMaj = 6\n",
    "dMin = 4\n",
    "p = 0.5\n",
    "\n",
    "\n",
    "# Calculation of q_LR (variable qbar in this code)\n",
    "\n",
    "qvec = range(0.501,length=1000,stop=0.999)\n",
    "z = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    z[k] = log((gammahat*(1-qvec[k])+(1-gammahat))/(gammahat*qvec[k]+(1-gammahat)))/log((1-qvec[k])/qvec[k])\n",
    "end\n",
    "qb = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    qb[k] = abs(1/2+((gamma-(2-gamma)*z[k])*(dMaj-dMin))/(2*gamma*(2*(1+n)+(dMaj+dMin)*(1+z[k])))-qvec[k])\n",
    "end\n",
    "diff = 0\n",
    "ind = 0\n",
    "diff, ind = findmin(qb)\n",
    "qbar = qvec[ind]\n",
    "\n",
    "\n",
    "# Determination of q\n",
    "\n",
    "q = 0.3*0.5+0.7*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beliefs with misperception\n",
    "\n",
    "n = 4+1\n",
    "dMaj = 6\n",
    "dMin = 4\n",
    "\n",
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "epsilon = 0.03\n",
    "\n",
    "for i in 1:1:4600\n",
    "    mu1vec_m[1,i] = mu1vec_m[1,i]+epsilon\n",
    "    mu2vec_m[1,i] = mu2vec_m[1,i]+epsilon\n",
    "end\n",
    "\n",
    "for i in 4601:1:10000\n",
    "    mu1vec_m[1,i] = mu1vec_m[1,i]-epsilon\n",
    "    mu2vec_m[1,i] = mu2vec_m[1,i]-epsilon\n",
    "end\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:50\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\global_warming_simulation.png\""
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share1vec_m = zeros(T+1)\n",
    "share1vec_m[1] = 0.46\n",
    "share2vec_m = zeros(T+1)\n",
    "share2vec_m[1] = 0.46\n",
    "\n",
    "for i in 1:1:10000\n",
    "    for t in 1:1:T\n",
    "        if mu1vec_m[t+1,i]>0.5\n",
    "            share1vec_m[t+1] = share1vec_m[t+1]+1/10000\n",
    "        end\n",
    "        if mu2vec_m[t+1,i]>0.5\n",
    "            share2vec_m[t+1] = share2vec_m[t+1]+1/10000\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "plot_data = [share1vec_m share2vec_m]\n",
    "plot_labels = [\"Democrats\" \"Republicans\"]\n",
    "plot_colors = [:blue :red]\n",
    "plot_styles = [:solid :solid]\n",
    "plot(plot_data, label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), dpi=600, legend=:topright, ylabel = \"Percent of population\", xlabel = \"Period, t\", legendfontsize=12, labelfontsize=14, tickfontsize=12)\n",
    "savefig(\"global_warming_simulation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5041894894894895"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simulation of climate change opinions (Gallup) with expanding friend networks\n",
    "\n",
    "T = 50\n",
    "n = 5\n",
    "dMaj = 6\n",
    "dMin = 4\n",
    "p = 0.5\n",
    "\n",
    "\n",
    "# Calculation of q_LR (variable qbar in this code)\n",
    "\n",
    "qvec = range(0.501,length=1000,stop=0.999)\n",
    "z = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    z[k] = log((gammahat*(1-qvec[k])+(1-gammahat))/(gammahat*qvec[k]+(1-gammahat)))/log((1-qvec[k])/qvec[k])\n",
    "end\n",
    "qb = zeros(100,1)\n",
    "for k in 1:1:100\n",
    "    qb[k] = abs(1/2+((gamma-(2-gamma)*z[k])*(dMaj-dMin))/(2*gamma*(2*(1+n)+(dMaj+dMin)*(1+z[k])))-qvec[k])\n",
    "end\n",
    "diff = 0\n",
    "ind = 0\n",
    "diff, ind = findmin(qb)\n",
    "qbar = qvec[ind]\n",
    "\n",
    "\n",
    "# Determination of q\n",
    "\n",
    "q = 0.3*0.5+0.7*qbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beliefs with misperception\n",
    "\n",
    "n = 4+1\n",
    "dMaj = 6\n",
    "dMin = 4\n",
    "\n",
    "mu1vec_m = zeros(T+1,10000)\n",
    "mu1vec_m[1,:] = p*ones(1,10000)\n",
    "mu2vec_m = zeros(T+1,10000)\n",
    "mu2vec_m[1,:] = p*ones(1,10000)\n",
    "\n",
    "epsilon = 0.03\n",
    "\n",
    "for i in 1:1:4600\n",
    "    mu1vec_m[1,i] = mu1vec_m[1,i]+epsilon\n",
    "    mu2vec_m[1,i] = mu2vec_m[1,i]+epsilon\n",
    "end\n",
    "\n",
    "for i in 4601:1:10000\n",
    "    mu1vec_m[1,i] = mu1vec_m[1,i]-epsilon\n",
    "    mu2vec_m[1,i] = mu2vec_m[1,i]-epsilon\n",
    "end\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 1:1:25\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end\n",
    "\n",
    "n = 4*4+1\n",
    "dMaj = 6*4\n",
    "dMin = 4*4\n",
    "\n",
    "for i in 1:1:10000\n",
    "for t in 26:1:50\n",
    "    sMaj1 = rand(Binomial(dMaj,gamma))\n",
    "    sMin1 = rand(Binomial(dMin,gamma))\n",
    "    sn1 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj1 = rand(Binomial(sMaj1,q))\n",
    "    sAmin1 = rand(Binomial(sMin1,q))\n",
    "    sAn1 = rand(Binomial(sn1,q))\n",
    "    \n",
    "    mu1vec_m[t+1,i] = mu1vec_m[t,i]/(mu1vec_m[t,i]+(1-mu1vec_m[t,i])*((1-q)/q)^(sAmaj1+sAn1-(sMin1-sAmin1+sn1-sAn1))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMin-sMin1+sAmin1)-(dMaj-sAmaj1)))\n",
    "    \n",
    "    sMaj2 = rand(Binomial(dMin,gamma))\n",
    "    sMin2 = rand(Binomial(dMaj,gamma))\n",
    "    sn2 = rand(Binomial(n,gamma))\n",
    "    \n",
    "    sAmaj2 = rand(Binomial(sMaj2,q))\n",
    "    sAmin2 = rand(Binomial(sMin2,q))\n",
    "    sAn2 = rand(Binomial(sn2,q))\n",
    "    \n",
    "    mu2vec_m[t+1,i] = mu2vec_m[t,i]/(mu2vec_m[t,i]+(1-mu2vec_m[t,i])*((1-q)/q)^(sAmaj2+sAn2-(sMin2-sAmin2+sn2-sAn2))*((gammahat*(1-q)+(1-gammahat))/(gammahat*q+(1-gammahat)))^((dMaj-sMin2+sAmin2)-(dMin-sAmaj2)))\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"C:\\\\Users\\\\danil\\\\Documents\\\\Julia code\\\\global_warming_simulation_expansion.png\""
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share1vec_m = zeros(T+1)\n",
    "share1vec_m[1] = 0.46\n",
    "share2vec_m = zeros(T+1)\n",
    "share2vec_m[1] = 0.46\n",
    "\n",
    "for i in 1:1:10000\n",
    "    for t in 1:1:T\n",
    "        if mu1vec_m[t+1,i]>0.5\n",
    "            share1vec_m[t+1] = share1vec_m[t+1]+1/10000\n",
    "        end\n",
    "        if mu2vec_m[t+1,i]>0.5\n",
    "            share2vec_m[t+1] = share2vec_m[t+1]+1/10000\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "plot_data = [share1vec_m share2vec_m]\n",
    "plot_labels = [\"Democrats\" \"Republicans\"]\n",
    "plot_colors = [:blue :red]\n",
    "plot_styles = [:solid :solid]\n",
    "plot(plot_data, label = plot_labels, linecolor = plot_colors, linestyle = plot_styles, ylims=[0,1], size=(800,600), dpi=600, legend=:topright, ylabel = \"Percent of population\", xlabel = \"Period, t\", legendfontsize=12, labelfontsize=14, tickfontsize=12)\n",
    "savefig(\"global_warming_simulation_expansion\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
